<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>cgglse.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CGGLSE.1"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</a>( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      COMPLEX            A( LDA, * ), B( LDB, * ), C( * ), D( * ),
     $                   WORK( * ), X( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="CGGLSE.19"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</a> solves the linear equality-constrained least squares (LSE)
</span><span class="comment">*</span><span class="comment">  problem:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          minimize || c - A*x ||_2   subject to   B*x = d
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
</span><span class="comment">*</span><span class="comment">  M-vector, and d is a given P-vector. It is assumed that
</span><span class="comment">*</span><span class="comment">  P &lt;= N &lt;= M+P, and
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           rank(B) = P and  rank( (A) ) = N.
</span><span class="comment">*</span><span class="comment">                                ( (B) )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  These conditions ensure that the LSE problem has a unique solution,
</span><span class="comment">*</span><span class="comment">  which is obtained using a generalized RQ factorization of the
</span><span class="comment">*</span><span class="comment">  matrices (B, A) given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     B = (0 R)*Q,   A = Z*T*Q.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of rows of the matrix A.  M &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrices A and B. N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  P       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of rows of the matrix B. 0 &lt;= P &lt;= N &lt;= M+P.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) COMPLEX array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the M-by-N matrix A.
</span><span class="comment">*</span><span class="comment">          On exit, the elements on and above the diagonal of the array
</span><span class="comment">*</span><span class="comment">          contain the min(M,N)-by-N upper trapezoidal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A. LDA &gt;= max(1,M).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) COMPLEX array, dimension (LDB,N)
</span><span class="comment">*</span><span class="comment">          On entry, the P-by-N matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
</span><span class="comment">*</span><span class="comment">          contains the P-by-P upper triangular matrix R.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B. LDB &gt;= max(1,P).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input/output) COMPLEX array, dimension (M)
</span><span class="comment">*</span><span class="comment">          On entry, C contains the right hand side vector for the
</span><span class="comment">*</span><span class="comment">          least squares part of the LSE problem.
</span><span class="comment">*</span><span class="comment">          On exit, the residual sum of squares for the solution
</span><span class="comment">*</span><span class="comment">          is given by the sum of squares of elements N-P+1 to M of
</span><span class="comment">*</span><span class="comment">          vector C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) COMPLEX array, dimension (P)
</span><span class="comment">*</span><span class="comment">          On entry, D contains the right hand side vector for the
</span><span class="comment">*</span><span class="comment">          constrained equation.
</span><span class="comment">*</span><span class="comment">          On exit, D is destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  X       (output) COMPLEX array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On exit, X is the solution of the LSE problem.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK. LWORK &gt;= max(1,M+N+P).
</span><span class="comment">*</span><span class="comment">          For optimum performance LWORK &gt;= P+min(M,N)+max(M,N)*NB,
</span><span class="comment">*</span><span class="comment">          where NB is an upper bound for the optimal blocksizes for
</span><span class="comment">*</span><span class="comment">          <a name="CGEQRF.87"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>, <a name="CGERQF.87"></a><a href="cgerqf.f.html#CGERQF.1">CGERQF</a>, <a name="CUNMQR.87"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a> and <a name="CUNMRQ.87"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.92"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">          = 1:  the upper triangular factor R associated with B in the
</span><span class="comment">*</span><span class="comment">                generalized RQ factorization of the pair (B, A) is
</span><span class="comment">*</span><span class="comment">                singular, so that rank(B) &lt; P; the least squares
</span><span class="comment">*</span><span class="comment">                solution could not be computed.
</span><span class="comment">*</span><span class="comment">          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
</span><span class="comment">*</span><span class="comment">                T associated with A in the generalized RQ factorization
</span><span class="comment">*</span><span class="comment">                of the pair (B, A) is singular, so that
</span><span class="comment">*</span><span class="comment">                rank( (A) ) &lt; N; the least squares solution could not
</span><span class="comment">*</span><span class="comment">                    ( (B) )
</span><span class="comment">*</span><span class="comment">                be computed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      COMPLEX            CONE
      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY
      INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
     $                   NB4, NR
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CAXPY, CCOPY, CGEMV, <a name="CGGRQF.120"></a><a href="cggrqf.f.html#CGGRQF.1">CGGRQF</a>, CTRMV, <a name="CTRTRS.120"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>,
     $                   <a name="CUNMQR.121"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>, <a name="CUNMRQ.121"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>, <a name="XERBLA.121"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.124"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      EXTERNAL           <a name="ILAENV.125"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> 
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          INT, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      MN = MIN( M, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Calculate workspace
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0) THEN
         IF( N.EQ.0 ) THEN
            LWKMIN = 1
            LWKOPT = 1
         ELSE
            NB1 = <a name="ILAENV.156"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="CGEQRF.156"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            NB2 = <a name="ILAENV.157"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="CGERQF.157"></a><a href="cgerqf.f.html#CGERQF.1">CGERQF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            NB3 = <a name="ILAENV.158"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="CUNMQR.158"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>'</span>, <span class="string">' '</span>, M, N, P, -1 )
            NB4 = <a name="ILAENV.159"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="CUNMRQ.159"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>'</span>, <span class="string">' '</span>, M, N, P, -1 )
            NB = MAX( NB1, NB2, NB3, NB4 )
            LWKMIN = M + N + P
            LWKOPT = P + MN + MAX( M, N )*NB
         END IF
         WORK( 1 ) = LWKOPT
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
            INFO = -12
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.172"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGGLSE.172"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the GRQ factorization of matrices B and A:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            B*Q' = (  0  T12 ) P   Z'*A*Q' = ( R11 R12 ) N-P
</span><span class="comment">*</span><span class="comment">                     N-P  P                  (  0  R22 ) M+P-N
</span><span class="comment">*</span><span class="comment">                                               N-P  P
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     where T12 and R11 are upper triangular, and Q and Z are
</span><span class="comment">*</span><span class="comment">     unitary.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGGRQF.192"></a><a href="cggrqf.f.html#CGGRQF.1">CGGRQF</a>( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
     $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
      LOPT = WORK( P+MN+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update c = Z'*c = ( c1 ) N-P
</span><span class="comment">*</span><span class="comment">                       ( c2 ) M+P-N
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMQR.199"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate Transpose'</span>, M, 1, MN, A, LDA,
     $             WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ),
     $             LWORK-P-MN, INFO )
      LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve T12*x2 = d for x2
</span><span class="comment">*</span><span class="comment">
</span>      IF( P.GT.0 ) THEN
         CALL <a name="CTRTRS.207"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, P, 1,
     $                B( 1, N-P+1 ), LDB, D, P, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 1
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Put the solution in X
</span><span class="comment">*</span><span class="comment">
</span>      CALL CCOPY( P, D, 1, X( N-P+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update c1
</span><span class="comment">*</span><span class="comment">
</span>         CALL CGEMV( <span class="string">'No transpose'</span>, N-P, P, -CONE, A( 1, N-P+1 ), LDA,
     $               D, 1, CONE, C, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve R11*x1 = c1 for x1
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.GT.P ) THEN
         CALL <a name="CTRTRS.228"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non-unit'</span>, N-P, 1,
     $                A, LDA, C, N-P, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 2
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Put the solutions in X
</span><span class="comment">*</span><span class="comment">
</span>         CALL CCOPY( N-P, C, 1, X, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the residual vector:
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LT.N ) THEN
         NR = M + P - N
         IF( NR.GT.0 )
     $      CALL CGEMV( <span class="string">'No transpose'</span>, NR, N-M, -CONE, A( N-P+1, M+1 ),
     $                  LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 )
      ELSE
         NR = P
      END IF
      IF( NR.GT.0 ) THEN
         CALL CTRMV( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non unit'</span>, NR,
     $               A( N-P+1, N-P+1 ), LDA, D, 1 )
         CALL CAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Backward transformation x = Q'*x
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMRQ.259"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate Transpose'</span>, N, 1, P, B, LDB,
     $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
      WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CGGLSE.265"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
